1 Overview

This project aims to connect floral phenology to reproductive success on the individual level for three subalpine plants. Data were collected at the Rocky Mountain Biological Laboratory during the 2019 field season. Focal species include Mertensia fusiformis, Potentilla pulcherrima, and Delphinium nuttallianum.

Our questions specifically ask: 1. how does blooming time affect seed set for individiuals in a population? 2. how do abiotic factors such as soil moisture affect individual seed set? 3. how do biotic factors such as conspecific density and pollen limitation affect individual seed set?

We tagged individuals of the three species in control and manipulated plots in eight sites across East River valley and Washington Gulch. We tracked blooming time of those tagged individuals and conducted a pollen limitation experiment with two treatments, open and pollen-supplemented. Seed set was quantified as the count of total developed seeds and proportion of developed seeds in each individual.

Because there were sources of non-independence in this study, our statistical approach involves generalized linear mixed effects models. Our models includes site as a random effect and the following fixed effects: the interaction between being before or after the population peak and deviation from the peak and the interaction between conspecific density and plot treatment. A separate model included pollen limitation treatment.

2 Load Packages

rm(list = ls()) #clearing environment

suppressWarnings(suppressPackageStartupMessages(require(knitr)))
suppressWarnings(suppressPackageStartupMessages(require(tidyverse)))
suppressWarnings(suppressPackageStartupMessages(require(lme4)))
suppressWarnings(suppressPackageStartupMessages(require(ggplot2)))
suppressWarnings(suppressPackageStartupMessages(require(RColorBrewer)))
suppressWarnings(suppressPackageStartupMessages(require(fitdistrplus)))
suppressWarnings(suppressPackageStartupMessages(require(glmmADMB)))
suppressWarnings(suppressPackageStartupMessages(require(MuMIn)))
suppressWarnings(suppressPackageStartupMessages(require(lubridate)))
suppressWarnings(suppressPackageStartupMessages(require(DHARMa)))
suppressWarnings(suppressPackageStartupMessages(require(glmmTMB)))
suppressWarnings(suppressPackageStartupMessages(require(dplyr)))
suppressWarnings(suppressPackageStartupMessages(require(stringr)))
suppressWarnings(suppressPackageStartupMessages(require(scales)))

3 Data Import & Basic Formatting

The imported csv files contain the following for each individual of a species: bloom start and end estimates, relative position of blooming in the community, soil moisture variables, and seed counts. These elements were compiled from five different data sets and calculated in a separate document. To view the data cleaning process, refer to the data cleaning R Markdown file, Data_Cleaning.Rmd.

3.1 Mertensia Data

The first dataset contains the individual data for Mertensia fusiformis, which was found in four of our sites. For Mertensia, we calculated the proportion of developed seeds in addition to the seed counts. We also reformatted the relative position as a deviation variable (the absolute value of the relative position) and a early/late variable (the direction of the relative position).

mert = read.csv("Mertensia_final2.csv") # read in Mertensia data
mert<-mert[!is.na(mert$relative.position),] #removing NA values from relative position column
mert<-mert[!is.na(mert$dev.seed),]
mert$prop.dev.seeds<-(mert$dev.seed/mert$total.seeds) #adding proportion of developed seeds
mert$deviation<-abs(mert$relative.position) #getting deviation from the peak
mert$early.late<-factor(NA, levels = c("early","late")) #creating early/late variables based on relative position
mert$early.late[mert$relative.position<0]<-"early"
mert$early.late[mert$relative.position>0]<-"late"
mert<-mert[!is.na(mert$early.late),]

3.2 Delphinium Data

The second dataset contains individual data for Delphinium nuttallianum. Delphinium was found in four of our sites. For Delphinium, we calculated proportion, deviation, and early/late value like we did with the Mertensia data.

delph = read.csv("Delphinium_final2.csv") # read in Delphinium data
delph<-delph[!is.na(delph$total.seeds),] #removing NA values in seed column
delph<-delph[!is.na(delph$relative.position),] #removing NA values in the relative position
delph$undev.seed1[is.na(delph$undev.seed1)] <- 0 #converting NA to 0
delph$prop.dev.seeds<-(delph$total.seeds/(delph$undev.seed1+delph$undev.seed+delph$total.seeds)) #adding proportion of developed seeds
delph$deviation<-abs(delph$relative.position) #getting deviation from the peak
delph$early.late<-factor(NA, levels = c("early","late")) #creating early/late variables based on relative position
delph$early.late[delph$relative.position<0]<-"early"
delph$early.late[delph$relative.position>0]<-"late"
delph<-delph[!is.na(delph$early.late),]

3.3 Potentilla Data

The third dataset contains Potentilla pulcherrima data. Potentilla was set up in all but one of our sites. Potentilla data only includes total seed counts because we couldn’t count the undeveloped seeds. However, we calculated the deviation and early/late values like the previous species.

pot = read.csv("Potentilla_final2.csv") # read in Potentilla data
pot<-pot[c(-119,-223),] #removing one row with NA in treatment column
pot<-pot[!is.na(pot$relative.position),] #removing NA values in relative position
pot<-pot[!is.na(pot$total.seeds),] #removing NA values in total seeds column
pot$deviation<-abs(pot$relative.position) #getting deviation from the peak
pot$early.late<-factor(NA, levels = c("early","late")) #creating early/late variables based on relative position
pot$early.late[pot$relative.position<0]<-"early"
pot$early.late[pot$relative.position>0]<-"late"
pot<-pot[!is.na(pot$early.late),]

Because many of the Potentilla individuals were still in bloom when we finished our field season, we had to calculate a midpoint and relative position based on those individuals that had senesced. A key factor that determines duration of blooming is the number of flowers on the individual. The data imported includes those relative position estimates for the individuals that didn’t senesce before we finished our field season. For more information on how those numbers were calculated, see Data_Cleaning.Rmd.

3.4 Flower Count Data

The community flower count data is essential to the number of conspecifics effect. Each of the species has to be analyzed separately.

flower.counts=read.csv("flower_counts.csv") #read in flower count data

3.4.1 Mertensia

We calculated the average number of Mertensia flowers in each site and plot in each week. This number was then merged with the Mertensia individuals to add the number of conspecifics variable.

mert.flower.counts<-flower.counts[(flower.counts$plant=="Mertensia fusiformis"),] #limit flower counts to Mertensia data
mert.counts.by.site<-aggregate(x=mert.flower.counts$total_flowers,by=list(mert.flower.counts$date,mert.flower.counts$site,mert.flower.counts$plot),FUN="mean") #calculate mean flower counts by week and site
colnames(mert.counts.by.site)<-c("date","site","plot.treat","number.conspecifics")

mert.counts.by.site$site<-as.character(mert.counts.by.site$site) #changing class to adjust names
mert.counts.by.site$plot.treat<-as.character(mert.counts.by.site$plot.treat)

mert.counts.by.site$site[mert.counts.by.site$site=="Avery Picnic"]<-"AveryPicnic" #fixing names
mert.counts.by.site$site[mert.counts.by.site$site=="Stupid Falls"]<-"StupidFalls"
mert.counts.by.site$site[mert.counts.by.site$site=="Wash 3C"]<-"WG3C"
mert.counts.by.site$plot.treat[!(mert.counts.by.site$plot.treat=="Kaysee")]
mert.counts.by.site$plot.treat[mert.counts.by.site$plot.treat=="Control"]<-"control"
mert.counts.by.site$plot.treat[mert.counts.by.site$plot.treat=="Manipulated"]<-"manip"

mert$week.estimate<-strftime(mert$midpoint, format = "%V") #calculating week numbers in preparation for the merge
mert.counts.by.site$week.estimate<-strftime(mert.counts.by.site$date, format = "%V")

mert<-merge(mert,mert.counts.by.site, by.x = c("site","plot.treat","week.estimate"), by.y = c("site","plot.treat","week.estimate")) #merging by site, plot treatment, and week estimate

3.4.2 Delphinium

The same was repeated for Delphinium.

delph.flower.counts<-flower.counts[(flower.counts$plant=="Delphinium nuttallianum"),] #limit flower counts to Delphinium data
delph.counts.by.site<-aggregate(x=delph.flower.counts$total_flowers,by=list(delph.flower.counts$date,delph.flower.counts$site,delph.flower.counts$plot),FUN="mean") #calculate mean flower counts by week and site
colnames(delph.counts.by.site)<-c("date","site","plot.treat","number.conspecifics")

delph.counts.by.site$site<-as.character(delph.counts.by.site$site) #changing class to adjust names
delph.counts.by.site$plot.treat<-as.character(delph.counts.by.site$plot.treat)

delph.counts.by.site$site[delph.counts.by.site$site=="Rustlers Gulch"]<-"RustlersGulch" #fixing names
delph.counts.by.site$site[delph.counts.by.site$site=="Wash 3C"]<-"WG3C"
#delph.counts.by.site$plot.treat[!(delph.counts.by.site$plot.treat=="Kaysee")]
delph.counts.by.site$plot.treat[delph.counts.by.site$plot.treat=="Control"]<-"control"
delph.counts.by.site$plot.treat[delph.counts.by.site$plot.treat=="Manipulated"]<-"manip"

delph$week.estimate<-strftime(delph$midpoint, format = "%V") #calculating week numbers in preparation for the merge
delph.counts.by.site$week.estimate<-strftime(delph.counts.by.site$date, format = "%V")

delph<-merge(delph,delph.counts.by.site, by.x = c("site","plot.treat","week.estimate"), by.y = c("site","plot.treat","week.estimate")) #merging by site, plot treatment, and week estimate

3.4.3 Potentilla

The same was done for Potentilla.

pot.flower.counts<-flower.counts[(flower.counts$plant=="Potentilla pulcherrima"),] #limit flower counts to Potentilla data
pot.counts.by.site<-aggregate(x=pot.flower.counts$total_flowers,by=list(pot.flower.counts$date,pot.flower.counts$site,pot.flower.counts$plot),FUN="mean") #calculate mean flower counts by week and site
colnames(pot.counts.by.site)<-c("date","site","plot.treat","number.conspecifics")

pot.counts.by.site$site<-as.character(pot.counts.by.site$site) #changing class to adjust names
pot.counts.by.site$plot.treat<-as.character(pot.counts.by.site$plot.treat)

pot.counts.by.site$site[pot.counts.by.site$site=="Rustlers Gulch"]<-"RustlersGulch" #fixing names
pot.counts.by.site$site[pot.counts.by.site$site=="Wash 3C"]<-"WG3C"
pot.counts.by.site$site[pot.counts.by.site$site=="Bellview Bench"]<-"BellviewBench"
pot.counts.by.site$site[pot.counts.by.site$site=="Avery Picnic"]<-"AveryPicnic"
pot.counts.by.site$site[pot.counts.by.site$site=="Stupid Falls"]<-"StupidFalls"
#pot.counts.by.site$plot.treat[!(pot.counts.by.site$plot.treat=="Kaysee")]
pot.counts.by.site$plot.treat[pot.counts.by.site$plot.treat=="Control"]<-"control"
pot.counts.by.site$plot.treat[pot.counts.by.site$plot.treat=="Manipulated"]<-"manip"

pot$week.estimate<-strftime(pot$midpoint, format = "%V") #calculating week numbers in preparation for the merge
pot.counts.by.site$week.estimate<-strftime(pot.counts.by.site$date, format = "%V")

pot<-merge(pot,pot.counts.by.site, by.x = c("site","plot.treat","week.estimate"), by.y = c("site","plot.treat","week.estimate")) #merging by site, plot treatment, and week estimate

4 Mixed Effects Models: Effect of Soil Moisture Variables on Total Seeds

Soil moisture is a covariate that may impact seed production. In the field, we measured soil moisture in designated areas of the plot once per week, for 9 weeks. To understand which variables representing change in soil moisture are most influential for seed production, we ran a generalized linear mixed model with all of the soil moisture variables. Site, plot treatment, and species are random effects, while the soil moisture variables (mean moisture, rate of change, effective minimum or value at week 7, final moisture reading, and range) are fixed effects. Once we had this model, we ran the dredge function to find the most appropriate combination of variables for our larger model.

We are including all of the individuals of the three species and using species as a random effect.

condensed.mert<-mert[,c(-4,-20,-21,-23,-24,-30)] #removing unnecessary columns
condensed.delph<-delph[,c(-4,-20:-28,-30,-36)]
condensed.pot<-pot[,c(-4,-20:-28,-30,-33:-37)]
names(condensed.mert)[names(condensed.mert)=="dev.seed"]<-"total.seeds"

condensed.mert<-condensed.mert[(condensed.mert$treat=="open"),] #limiting to open treatment
condensed.delph<-condensed.delph[(condensed.delph$treat=="open"),] 
condensed.pot<-condensed.pot[(condensed.pot$treat=="open"),] 

all.seeds<-rbind(condensed.mert,condensed.delph,condensed.pot) #combining all of the data sets

all.seed.nbinom<-fitdist(all.seeds$total.seeds, "nbinom") #fitting Mertensia data to negative binomial distribution using fitdistplus
plot(all.seed.nbinom) #plotting to see the fit

The seeds fit a negative binomial distribution, so we used the glmmTMB package.

soil.mod<-glmmTMB(total.seeds ~ mean + rate + effect.min + end + range + (1|site/plot.treat) + (1|species), family = "nbinom2", data = all.seeds) #model for soil moisture variables with all seeds
summary(soil.mod)

##  Family: nbinom2  ( log )
## Formula:          
## total.seeds ~ mean + rate + effect.min + end + range + (1 | site/plot.treat) +  
##     (1 | species)
## Data: all.seeds
## 
##      AIC      BIC   logLik deviance df.resid 
##   1238.8   1266.8   -609.4   1218.8      112 
## 
## Random effects:
## 
## Conditional model:
##  Groups          Name        Variance  Std.Dev. 
##  plot.treat:site (Intercept) 8.441e-08 0.0002905
##  site            (Intercept) 4.806e-02 0.2192257
##  species         (Intercept) 5.434e-01 0.7371681
## Number of obs: 122, groups:  plot.treat:site, 14; site, 7; species, 3
## 
## Overdispersion parameter for nbinom2 family (): 1.34 
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.995718   0.664878   4.506 6.62e-06 ***
## mean        -0.017207   0.083280  -0.207    0.836    
## rate        -0.007747   0.159466  -0.049    0.961    
## effect.min   0.082413   0.101464   0.812    0.417    
## end          0.030435   0.042880   0.710    0.478    
## range        0.029783   0.083084   0.358    0.720    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

dredge_result<-dredge(soil.mod) #using dredge function for model selection
subset(dredge_result, delta<4)

## Global model call: glmmTMB(formula = total.seeds ~ mean + rate + effect.min + end + 
##     range + (1 | site/plot.treat) + (1 | species), data = all.seeds, 
##     family = "nbinom2", ziformula = ~0, dispformula = ~1)
## ---
## Model selection table 
##    cnd((Int)) dsp((Int)) cnd(eff.min) cnd(end) cnd(men)  cnd(rng) cnd(rat)
## 17      3.886          +                                          -0.07427
## 1       4.015          +                                                  
## 18      3.323          +      0.06957                             -0.07516
## 10      3.101          +      0.07761                    0.029490         
## 3       3.428          +               0.05658                            
## 2       3.475          +      0.06755                                     
## 9       3.765          +                                 0.025740         
## 5       3.456          +                        0.03591                   
## 21      3.416          +                        0.03126           -0.06544
## 19      3.455          +               0.04396                    -0.05917
## 4       3.114          +      0.05359  0.04546                            
## 11      3.329          +               0.04691           0.020350         
## 7       3.183          +               0.04490  0.02351                   
## 20      3.121          +      0.05915  0.02933                    -0.06408
## 12      2.929          +      0.06616  0.02979           0.024740         
## 26      3.183          +      0.07453                    0.015930 -0.04322
## 22      3.195          +      0.05685           0.01524           -0.07120
## 6       3.299          +      0.05036           0.02000                   
## 14      3.167          +      0.09727          -0.01950  0.037530         
## 25      3.844          +                                 0.006873 -0.06031
## 13      3.464          +                        0.02408  0.018220         
## 29      3.307          +                        0.05257 -0.034310 -0.12990
## 23      3.222          +               0.03227  0.02309           -0.05682
## 27      3.421          +               0.04371           0.005754 -0.04792
## 8       3.041          +      0.04554  0.04188  0.01113                   
## 15      3.209          +               0.04221  0.01314  0.016700         
## 16      2.987          +      0.08640  0.03120 -0.02082  0.033650         
## 28      3.016          +      0.06399  0.02736           0.013630 -0.03696
##    df   logLik   AICc delta weight
## 17  6 -610.872 1234.5  0.00  0.071
## 1   5 -611.994 1234.5  0.03  0.070
## 18  7 -609.782 1234.5  0.07  0.068
## 10  7 -609.799 1234.6  0.11  0.067
## 3   6 -610.938 1234.6  0.13  0.066
## 2   6 -611.002 1234.7  0.26  0.062
## 9   6 -611.141 1235.0  0.54  0.054
## 5   6 -611.207 1235.1  0.67  0.051
## 21  7 -610.196 1235.4  0.90  0.045
## 19  7 -610.241 1235.5  0.99  0.043
## 4   7 -610.350 1235.7  1.21  0.039
## 11  7 -610.423 1235.8  1.35  0.036
## 7   7 -610.630 1236.2  1.77  0.029
## 20  8 -609.519 1236.3  1.84  0.028
## 12  8 -609.527 1236.3  1.85  0.028
## 26  8 -609.628 1236.5  2.06  0.025
## 22  8 -609.638 1236.6  2.08  0.025
## 6   7 -610.800 1236.6  2.11  0.025
## 14  8 -609.682 1236.6  2.16  0.024
## 25  7 -610.847 1236.7  2.20  0.024
## 13  7 -610.872 1236.7  2.25  0.023
## 29  8 -609.869 1237.0  2.54  0.020
## 23  8 -609.895 1237.1  2.59  0.019
## 27  8 -610.222 1237.7  3.24  0.014
## 8   8 -610.286 1237.8  3.37  0.013
## 15  8 -610.347 1238.0  3.50  0.012
## 16  9 -609.379 1238.4  3.89  0.010
## 28  9 -609.399 1238.4  3.93  0.010
## Models ranked by AICc(x) 
## Random terms (all models): 
## 'cond(1 | site/plot.treat)', 'cond(1 | species)'

After running the model selection, we found that rate of soil moisture was the best model, but the null model was a very close second. The delta values between the rate of change predictor and the null model are 0.03 apart, which is very little. Even though the rate variable is not very different from the null model, we added the rate of change in our main models to incorporate the soil moisture covariate.

5 Mixed Effects Models with Count Data

Generalized linear mixed effects models were appropriate for this analysis because the data contain sources of non-independence. Individual plants were located in the same site as other individuals and were likely genetically similar. Therefore, we considered site a random effect.

The response variable is total seeds, which is discrete and non-zero. The interaction between deviation and the early/late factor was a fixed effect, as well as the interaction between plot treatment and conspecific density. We would have included soil moisture as a covariate and fixed effect in our model, but soil moisture was not important for seed production. Each species is analyzed separately.

5.1 Mertensia

We started by limiting the individuals to the open treatment (excluding the hand-pollinated treatment) and by checking the distribution of the total seed counts.

mert.open<-mert[(mert$treat=="open"),] #removed hand pollination treatment

mert.nbinom<-fitdist(mert.open$dev.seed, "nbinom") #fitting Mertensia data to negative binomial distribution using fitdistplus
plot(mert.nbinom) #plotting to see the fit

The Mertensia seed data fits a negative binomial distribution. The glmmTMB package is a good package for modeling data with a negative binomial. The glmmTMB package is also suitable for data that could be overdispersed or zero-inflated.

#mert.glmmtmb<-glmmTMB(dev.seed ~ deviation * early.late + plot.treat * number.conspecifics + rate + (1|site), family = "nbinom2", data = mert.open) #putting count data into glmmTMB for DHARMa package
#summary(mert.glmmtmb) #summary of glmmTMB model

We did not detect an effect of relative blooming time, conspecific density, or plot treatment on the total number of seeds in Mertensia individuals.

We then checked for overdispersion and zero-inflation in the data using the DHARMa package. The overdispersion and zero-inflation tests work with the glmmTMB package.

#mert.counts.simulation<-simulateResiduals(fittedModel = mert.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
#plot(mert.counts.simulation) #plotting simulation
#testDispersion(mert.counts.simulation) #checking for overdispersion
#testZeroInflation(mert.counts.simulation) #checking for zero inflation

The Mertensia count data is not overdispersed, but it is slightly zero-inflated. To address this, we adjusted our glmmtmb model to include zero-inflation.

mert.glmmtmb<-glmmTMB(dev.seed ~ deviation * early.late + plot.treat * number.conspecifics + rate + (1|site), ziformula=~1, family = "nbinom2", data = mert.open) #putting count data into glmmTMB for DHARMa package
summary(mert.glmmtmb) #summary of glmmTMB model

##  Family: nbinom2  ( log )
## Formula:          
## dev.seed ~ deviation * early.late + plot.treat * number.conspecifics +  
##     rate + (1 | site)
## Zero inflation:            ~1
## Data: mert.open
## 
##      AIC      BIC   logLik deviance df.resid 
##    397.5    417.3   -187.7    375.5       34 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.1723   0.4151  
## Number of obs: 45, groups:  site, 4
## 
## Overdispersion parameter for nbinom2 family (): 2.62 
## 
## Conditional model:
##                                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)                          3.633737   0.635848   5.715  1.1e-08
## deviation                            0.062103   0.092385   0.672    0.501
## early.latelate                       0.338458   0.503457   0.672    0.501
## plot.treatmanip                     -0.808112   0.613243  -1.318    0.188
## number.conspecifics                 -0.005894   0.004515  -1.305    0.192
## rate                                -0.055190   0.065045  -0.848    0.396
## deviation:early.latelate            -0.206822   0.172794  -1.197    0.231
## plot.treatmanip:number.conspecifics  0.012496   0.015982   0.782    0.434
##                                        
## (Intercept)                         ***
## deviation                              
## early.latelate                         
## plot.treatmanip                        
## number.conspecifics                    
## rate                                   
## deviation:early.latelate               
## plot.treatmanip:number.conspecifics    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -2.367      0.544  -4.351 1.36e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We then checked this new model for overdispersion and zero-inflation.

mert.counts.simulation<-simulateResiduals(fittedModel = mert.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
plot(mert.counts.simulation) #plotting simulation

testDispersion(mert.counts.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.92998, p-value = 0.968
## alternative hypothesis: two.sided

testZeroInflation(mert.counts.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0163, p-value = 1
## alternative hypothesis: two.sided

Zero-inflation is accounted for in the new model, and the model for Mertensia count data is no longer zero-inflated.

To test the pollen limitation individuals, we created a separate fixed effects model. This model has site and plot treatment as random effects and pollen limitation treatment as a fixed effect. The response variable is still total number of developed seeds in Mertensia individuals.

#mert.treat.glmmtmb<-glmmTMB(dev.seed ~ treat + (1|site/plot.treat), family = "nbinom2", data = mert) #Mertensia model for negative binomial data with glmmTMB package
#summary(mert.treat.glmmtmb) #summary of model output

We did not detect an effect of pollen limitation treatment on Mertensia total developed seeds.

We had to test the model assumptions of our pollen limitation model as well. After converting our model to an glmmTMB model, we tested for overdispersion and zero-inflation with the DHARMa package.

#mert.counts.treat.simulation<-simulateResiduals(fittedModel = mert.treat.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
#plot(mert.counts.treat.simulation) #plotting simulation
#testDispersion(mert.counts.treat.simulation) #checking for overdispersion
#testZeroInflation(mert.counts.treat.simulation) #checking for zero inflation

The data in the pollen limitation model for Mertensia is not overdispersed, but it is zero-inflated. We then adjusted our model to include zero-inflation and ran the overdispersion and zero-inflation tests again.

mert.treat.glmmtmb<-glmmTMB(dev.seed ~ treat + (1|site/plot.treat), ziformula=~1, family = "nbinom2", data = mert) #Mertensia model for negative binomial data with glmmTMB package
summary(mert.treat.glmmtmb) #summary of model output

##  Family: nbinom2  ( log )
## Formula:          dev.seed ~ treat + (1 | site/plot.treat)
## Zero inflation:            ~1
## Data: mert
## 
##      AIC      BIC   logLik deviance df.resid 
##    791.4    806.6   -389.7    779.4       87 
## 
## Random effects:
## 
## Conditional model:
##  Groups          Name        Variance Std.Dev.
##  plot.treat:site (Intercept) 0.08162  0.2857  
##  site            (Intercept) 0.01479  0.1216  
## Number of obs: 93, groups:  plot.treat:site, 8; site, 4
## 
## Overdispersion parameter for nbinom2 family (): 3.01 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    3.371413   0.155514  21.679   <2e-16 ***
## treatopen-hand 0.006098   0.133790   0.046    0.964    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.8927     0.4709  -6.143 8.12e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mert.counts.treat.simulation<-simulateResiduals(fittedModel = mert.treat.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
plot(mert.counts.treat.simulation) #plotting simulation

testDispersion(mert.counts.treat.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.96976, p-value = 0.984
## alternative hypothesis: two.sided

testZeroInflation(mert.counts.treat.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0187, p-value = 1
## alternative hypothesis: two.sided

The data are no longer zero-inflated, and we still did not detect an effect of pollen limitation on the Mertensia seed counts.

Below is the plot of seed set vs. relative position of blooming in the population for Mertensia individuals. Plot treatment is indicated by individual color.

ggplot(mert.open, aes(x=mert.open$relative.position, y=mert.open$dev.seed, group=mert.open$plot.treat)) +
  labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Mertensia Blooming Time on Total Developed Seeds") +
  geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
  scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
  geom_smooth()

We then plotted the interaction between deviation and early/late instead of relative position in order to fit a linear model.

ggplot(mert.open, aes(x=mert.open$deviation, y=mert.open$dev.seed, group= mert.open$early.late)) +
  labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Mertensia Deviation from Peak on Total Developed Seeds") + theme_classic() +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
  geom_smooth(method="glm",aes(color=early.late)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Before or After Peak") +
  scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))

Then we plotted the effect of conspecific density and plot treatment on total developed seeds.

ggplot(mert.open, aes(x=mert.open$number.conspecifics, y=mert.open$dev.seed, group= mert.open$plot.treat)) +
  labs(y="Total Developed Seeds",x="Conspecific Density", title="Effect of Mertensia Conspecific Density on Total Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
  geom_smooth(method="glm",aes(color=plot.treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))

Finally, we plotted all of the data with the pollinator active period. The pollinator active period was the time during which we saw interactions between pollinators and our species. The active period is shaded.

mert$midpoint2<-format(as.Date(mert$midpoint), format="%m-%d")

ggplot(mert, aes(x=mert$midpoint2, y=mert$dev.seed, group= mert$treat)) +
  labs(y="Total Developed Seeds",x="Date", title="Effect of Mertensia Blooming Time and Pollinator Activity on Total Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=treat)) +
  geom_smooth(aes(color=treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Open or Hand-pollinated") +
  scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  annotate("rect",xmin="06-19", xmax="07-03",ymin=0,ymax=Inf,alpha=0.1)

5.2 Delphinium

For Delphinium, we again limited the individuals to open treatment individuals.

delph.open<-delph[(delph$treat=="open"),] #removed hand pollination treatment

delph.nbinom<-fitdist(delph.open$total.seeds, "nbinom") #fitting negative binomial distribution using fitdistplus package
plot(delph.nbinom) #plotting distribution fit

The Delphinium data fits a negative binomial distribution. We used the glmmTMB package again.

delph.glmmtmb<-glmmTMB(total.seeds ~ deviation * early.late + plot.treat * number.conspecifics + rate + (1|site), family = "nbinom2", data = delph.open) #Delphinium model for negative binomial data with glmmTMB package
#removed treat fixed effect (only open)
summary(delph.glmmtmb) #summary of model output

##  Family: nbinom2  ( log )
## Formula:          
## total.seeds ~ deviation * early.late + plot.treat * number.conspecifics +  
##     rate + (1 | site)
## Data: delph.open
## 
##      AIC      BIC   logLik deviance df.resid 
##    409.5    426.6   -194.8    389.5       31 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  site   (Intercept) 8.014e-10 2.831e-05
## Number of obs: 41, groups:  site, 4
## 
## Overdispersion parameter for nbinom2 family (): 1.28 
## 
## Conditional model:
##                                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)                          2.999210   0.635426   4.720 2.36e-06
## deviation                            0.002673   0.115775   0.023    0.982
## early.latelate                       0.177613   0.822990   0.216    0.829
## plot.treatmanip                      0.952349   0.960652   0.991    0.322
## number.conspecifics                  0.001786   0.010991   0.162    0.871
## rate                                -0.141524   0.086406  -1.638    0.101
## deviation:early.latelate             0.161277   0.196370   0.821    0.411
## plot.treatmanip:number.conspecifics -0.009046   0.016771  -0.539    0.590
##                                        
## (Intercept)                         ***
## deviation                              
## early.latelate                         
## plot.treatmanip                        
## number.conspecifics                    
## rate                                   
## deviation:early.latelate               
## plot.treatmanip:number.conspecifics    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We did not detect an effect of relative blooming time, plot treatment, or conspecific density on the total developed seeds for Delphinium individuals.

Next, we tested the model assumptions as we did with Mertensia seed count data. We checked for overdispersion and zero-inflation.

delph.counts.simulation<-simulateResiduals(fittedModel = delph.glmmtmb) #creating simulation for Delphinium mixed effects model with DHARMa package
plot(delph.counts.simulation) #plotting simulation

testDispersion(delph.counts.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.83876, p-value = 0.536
## alternative hypothesis: two.sided

testZeroInflation(delph.counts.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 3.7594, p-value = 0.208
## alternative hypothesis: two.sided

The model for the Delphinium count data is not overdispersed or zero-inflated.

Again, we tested the pollen limitation in a separate mixed effects model.

#delph.treat.glmmtmb<-glmmTMB(total.seeds ~ treat + (1|site/plot.treat), family = "nbinom2", data = delph) #Mertensia model for negative binomial data with glmmTMB package
#summary(delph.treat.glmmtmb) #summary of model output

We did not detect an effect of pollen limitation on the total developed seeds of Delphinium individuals.

We again checked for overdispersion and zero-inflation in the pollen limitation data.

#delph.counts.treat.simulation<-simulateResiduals(fittedModel = delph.treat.glmmtmb) #creating simulation for Delphinium mixed effects model with DHARMa package
#plot(delph.counts.treat.simulation) #plotting simulation
#testDispersion(delph.counts.treat.simulation) #checking for overdispersion
#testZeroInflation(delph.counts.treat.simulation) #checking for zero inflation

The Delphinum pollen limitation model is not overdispersed. The data are slightly zero-inflated. Because the data were zero-inflated, we adjusted our model to include the zero-inflation formula and re-tested the overdispersion and zero-inflation.

delph.treat.glmmtmb<-glmmTMB(total.seeds ~ treat + (1|site/plot.treat), ziformula=~1, family = "nbinom2", data = delph) #Mertensia model for negative binomial data with glmmTMB package
summary(delph.treat.glmmtmb) #summary of model output

##  Family: nbinom2  ( log )
## Formula:          total.seeds ~ treat + (1 | site/plot.treat)
## Zero inflation:               ~1
## Data: delph
## 
##      AIC      BIC   logLik deviance df.resid 
##    824.1    838.7   -406.0    812.1       78 
## 
## Random effects:
## 
## Conditional model:
##  Groups          Name        Variance Std.Dev.
##  plot.treat:site (Intercept) 0.01460  0.1208  
##  site            (Intercept) 0.03166  0.1779  
## Number of obs: 84, groups:  plot.treat:site, 7; site, 4
## 
## Overdispersion parameter for nbinom2 family ():  1.6 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     3.87877    0.17023  22.786   <2e-16 ***
## treatopen-hand  0.09481    0.18614   0.509    0.611    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.0968     0.5684  -5.448 5.09e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

delph.counts.treat.simulation<-simulateResiduals(fittedModel = delph.treat.glmmtmb) #creating simulation for Delphinium mixed effects model with DHARMa package
plot(delph.counts.treat.simulation) #plotting simulation

testDispersion(delph.counts.treat.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.91298, p-value = 0.592
## alternative hypothesis: two.sided

testZeroInflation(delph.counts.treat.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0331, p-value = 1
## alternative hypothesis: two.sided

The new model is not overdispersed or zero-inflated. We still did not see an effect of pollen limitation on the Delphinium count data.

Below are the plots for the Delphinium individuals.

ggplot(delph.open, aes(x=delph.open$relative.position, y=delph.open$total.seeds, group=plot.treat)) +
  labs(title="Effect of Delphinium Blooming Time on Total Developed Seeds", y="Total Developed Seeds",x="Days Relative to Population Peak Bloom") +
  geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
  scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
  geom_smooth()

ggplot(delph.open, aes(x=delph.open$deviation, y=delph.open$total.seeds, group= delph.open$early.late)) +
  labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Delphinium Deviation from Peak on Total Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
  geom_smooth(method="glm",aes(color=early.late)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Before or After Peak") +
  scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))

ggplot(delph.open, aes(x=delph.open$number.conspecifics, y=delph.open$total.seeds, group= delph.open$plot.treat)) +
  labs(y="Total Developed Seeds",x="Conspecific Density", title="Effect of Delphinium Conspecific Density on Total Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
  geom_smooth(method="glm",aes(color=plot.treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))

Finally, we plotted the pollinator active period with the total number of developed seeds and midpoint of individuals.

delph$midpoint2<-format(as.Date(delph$midpoint), format="%m-%d")

ggplot(delph, aes(x=delph$midpoint2, y=delph$total.seeds, group= delph$treat)) +
  labs(y="Total Developed Seeds",x="Date", title="Effect of Delphinium Blooming Time and Pollinator Activity on Total Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=treat)) +
  geom_smooth(aes(color=treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Open or Hand-pollinated") +
  scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  annotate("rect",xmin="06-27", xmax="07-18",ymin=0,ymax=Inf,alpha=0.1)

5.3 Potentilla

We began by limiting the data to the open treatment Potentilla individuals.

pot.open<-pot[(pot$treat=="open"),] #removed hand pollination treatment

pot.nbinom<-fitdist(pot.open$total.seeds, "nbinom") #fitting negative binomial distribution using fitdistplus package
plot(pot.nbinom) #plotting distribution fit

The Potentilla seed data fits a negative binomial distribution. We used the glmmTMB package again.

#pot.glmmtmb<-glmmTMB(total.seeds ~ deviation * early.late + plot.treat * number.conspecifics + plot.treat + rate + (1|site), family = "nbinom2", data = pot.open) #Potentilla model for negative binomial data with glmmTMB package
#removed treat fixed effect (only open)
#summary(pot.glmmtmb) #summary of model output

We did not detect an effect of relative blooming time, plot treatment, or conspecific density on the total developed seeds of Potentilla individuals. However, the rate of change of soil moisture had an effect on seed set for Potentilla individuals. As rate of change increased, the total number of developed seeds decreased.

Below we checked overdispersion and zero-inflation.

#pot.counts.simulation<-simulateResiduals(fittedModel = pot.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
#plot(pot.counts.simulation) #plotting simulation
#testDispersion(pot.counts.simulation) #checking for overdispersion
#testZeroInflation(pot.counts.simulation) #checking for zero inflation

The Potentilla data are not overdispersed, but they are zero-inflated. To fix this, we had to include zero inflation in our model and again check for zero-inflation.

pot.glmmtmb<-glmmTMB(total.seeds ~ deviation * early.late + plot.treat * number.conspecifics + plot.treat + rate + (1|site), ziformula=~1, family = "nbinom2", data = pot.open) #new Potentilla model for negative binomial data with glmmTMB package
summary(pot.glmmtmb) #summary of model output

##  Family: nbinom2  ( log )
## Formula:          
## total.seeds ~ deviation * early.late + plot.treat * number.conspecifics +  
##     plot.treat + rate + (1 | site)
## Zero inflation:               ~1
## Data: pot.open
## 
##      AIC      BIC   logLik deviance df.resid 
##    433.8    451.2   -205.9    411.8       25 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  site   (Intercept) 7.935e-10 2.817e-05
## Number of obs: 36, groups:  site, 7
## 
## Overdispersion parameter for nbinom2 family (): 2.97 
## 
## Conditional model:
##                                       Estimate Std. Error z value Pr(>|z|)
## (Intercept)                          4.822e+00  3.842e-01  12.551  < 2e-16
## deviation                           -7.671e-03  2.079e-02  -0.369  0.71221
## early.latelate                      -7.089e-02  3.864e-01  -0.183  0.85444
## plot.treatmanip                      1.151e-01  3.615e-01   0.318  0.75019
## number.conspecifics                  5.871e-05  6.418e-03   0.009  0.99270
## rate                                -1.823e-01  6.961e-02  -2.619  0.00882
## deviation:early.latelate            -3.725e-02  3.157e-02  -1.180  0.23803
## plot.treatmanip:number.conspecifics -7.015e-04  6.932e-03  -0.101  0.91939
##                                        
## (Intercept)                         ***
## deviation                              
## early.latelate                         
## plot.treatmanip                        
## number.conspecifics                    
## rate                                ** 
## deviation:early.latelate               
## plot.treatmanip:number.conspecifics    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -3.557      1.016  -3.501 0.000463 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

When we incorporated zero-inflation into the model, we found a significant trend in the soil moisture variable. As the rate of change in soil moisture increases, the number of developed seeds in Potentilla individuals decreases. This was the same trend that we found in our previous model, but now the data are not overdispersed or zero-inflated.

However, we had to test for overdispersion and zero-inflation again.

pot.counts.simulation<-simulateResiduals(fittedModel = pot.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
plot(pot.counts.simulation) #plotting simulation

testDispersion(pot.counts.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.79932, p-value = 0.384
## alternative hypothesis: two.sided

testZeroInflation(pot.counts.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99206, p-value = 1
## alternative hypothesis: two.sided

The data are no longer overdispersed or zero-inflated for the Potentilla count data.

Below is the mixed effects model for the pollen limitation treatments.

#pot.treat.glmmtmb<-glmmTMB(total.seeds ~ treat + (1|site/plot.treat), family = "nbinom2", data = pot) #Potentilla model for negative binomial data with glmmadmb package
#summary(pot.treat.glmmtmb) #summary of model output

We did not detect an effect of the pollen limitation treatment on the total number of developed seeds for Potentilla individuals.

Below is the test for overdispersion and zero-inflation for Potentilla pollen limitation data.

#pot.counts.treat.simulation<-simulateResiduals(fittedModel = pot.treat.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
#plot(pot.counts.treat.simulation) #plotting simulation
#testDispersion(pot.counts.treat.simulation) #checking for overdispersion
#testZeroInflation(pot.counts.treat.simulation) #checking for zero inflation

The Potentilla pollen limitation data were not overdispersed, but they are zero-inflated. We addressed the zero-inflation issue by again including the zero-inflation formula into our models.

pot.treat.glmmtmb<-glmmTMB(total.seeds ~ treat + (1|site/plot.treat), ziformula=~1, family = "nbinom2", data = pot) #Potentilla model for negative binomial data with glmmadmb package
summary(pot.treat.glmmtmb) #summary of model output

##  Family: nbinom2  ( log )
## Formula:          total.seeds ~ treat + (1 | site/plot.treat)
## Zero inflation:               ~1
## Data: pot
## 
##      AIC      BIC   logLik deviance df.resid 
##    795.8    809.0   -391.9    783.8       61 
## 
## Random effects:
## 
## Conditional model:
##  Groups          Name        Variance  Std.Dev. 
##  plot.treat:site (Intercept) 1.274e-01 0.3569566
##  site            (Intercept) 2.947e-08 0.0001717
## Number of obs: 67, groups:  plot.treat:site, 11; site, 7
## 
## Overdispersion parameter for nbinom2 family (): 2.24 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     4.97899    0.16868  29.517   <2e-16 ***
## treatopen-hand -0.09409    0.17544  -0.536    0.592    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.0664     0.5943   -5.16 2.47e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pot.counts.treat.simulation<-simulateResiduals(fittedModel = pot.treat.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
plot(pot.counts.treat.simulation) #plotting simulation

testDispersion(pot.counts.treat.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 0.98533, p-value = 0.96
## alternative hypothesis: two.sided

testZeroInflation(pot.counts.treat.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.97656, p-value = 1
## alternative hypothesis: two.sided

The data are no longer zero-inflated, and we still did not detect an effect of pollen limitation in the Potentilla individuals.

Below are the plots for Potentilla individuals.

ggplot(pot.open, aes(x=pot.open$relative.position, y=pot.open$total.seeds, group=plot.treat)) +
  labs(title="Effect of Potentilla Blooming Time on Total Developed Seeds", y="Total Developed Seeds",x="Days Relative to Population Peak Bloom") +
  geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
  scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
  geom_smooth()

ggplot(pot.open, aes(x=pot.open$deviation, y=pot.open$total.seeds, group= pot.open$early.late)) +
  labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom",title="Effect of Potentilla Deviation from Peak on Total Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
  geom_smooth(method="glm",aes(color=early.late)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Before or After Peak") +
  scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))

ggplot(pot.open, aes(x=pot.open$number.conspecifics, y=pot.open$total.seeds, group= pot.open$plot.treat)) +
  labs(y="Total Developed Seeds",x="Conspecific Density", title="Effect of Potentilla Conspecific Density on Total Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
  geom_smooth(method="glm",aes(color=plot.treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))

Again we plotted the pollinator activity period.

ggplot(pot,aes(x=as.Date.factor(pot$midpoint), y=pot$total.seeds, group= pot$treat)) +
  labs(y="Total Developed Seeds",x="Date", title="Pollinator Activity Period and Total Developed Seeds in Potentilla") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=treat)) +
  geom_smooth(aes(color=treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Open or Hand-pollinated") +
  scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  annotate("rect",xmin=as.Date.character("2019-07-05"), xmax=as.Date.character("2019-08-24"),ymin=0,ymax=Inf,alpha=0.1) +
  scale_x_date(date_breaks = "3 days", date_labels = "%m/%d")

6 Mixed Effects Model with Proportion Data

Only the Mertensia and Delphinium data included proportion of developed seeds to total seeds. When counting seeds for the Potentilla individuals, undeveloped seeds were difficult to distinguish from other elements of the carpel and were not included. Because the data are proportions and we are using the cbind function, we assume a binomial distribution.

6.1 Mertensia

First, we had to calculate undeveloped seeds to use the cbind function for our proportion of developed seeds.

mert.open$undev.seeds<-(mert.open$total.seeds-mert.open$dev.seed) #creating an undeveloped seed column for cbind

Below is the mixed effects model for the proportion of developed seeds. The individuals are still only limited to the open treatment.

Mert.prop.glmmtmb<-glmmTMB(cbind(dev.seed,undev.seeds) ~ deviation * early.late + plot.treat * number.conspecifics + rate + (1|site), family = binomial, data = mert.open) #Mertensia model with proportion of developed seeds using the glmmTMB package
#removed treat fixed effect (only open in data set)
summary(Mert.prop.glmmtmb)

##  Family: binomial  ( logit )
## Formula:          
## cbind(dev.seed, undev.seeds) ~ deviation * early.late + plot.treat *  
##     number.conspecifics + rate + (1 | site)
## Data: mert.open
## 
##      AIC      BIC   logLik deviance df.resid 
##    299.2    315.4   -140.6    281.2       36 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  site   (Intercept) 0.03491  0.1868  
## Number of obs: 45, groups:  site, 4
## 
## Conditional model:
##                                      Estimate Std. Error z value Pr(>|z|)
## (Intercept)                          1.226241   0.302012   4.060  4.9e-05
## deviation                            0.005187   0.045449   0.114   0.9091
## early.latelate                      -0.100470   0.238530  -0.421   0.6736
## plot.treatmanip                     -0.421913   0.279412  -1.510   0.1310
## number.conspecifics                 -0.005196   0.002203  -2.358   0.0183
## rate                                -0.085739   0.035810  -2.394   0.0167
## deviation:early.latelate            -0.125525   0.082473  -1.522   0.1280
## plot.treatmanip:number.conspecifics -0.006530   0.007262  -0.899   0.3686
##                                        
## (Intercept)                         ***
## deviation                              
## early.latelate                         
## plot.treatmanip                        
## number.conspecifics                 *  
## rate                                *  
## deviation:early.latelate               
## plot.treatmanip:number.conspecifics    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A significant effect was detected in the number of conspecifics and the rate of change in soil moisture. As the number of conspecifics increased, the proportion of developed seeds in Mertensia individuals decreased. For the soil moisture variable, as the rate of change in soil moisture increased, the proportion of developed seeds decreased.

Once we had our model, we were able to check for overdispersion and zero inflaction. We checked for overdispersion and zero inflation with the DHARMa package.

mert.prop.simulation<-simulateResiduals(fittedModel = Mert.prop.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
plot(mert.prop.simulation) #plotting simulation

testDispersion(mert.prop.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.0336, p-value = 0.448
## alternative hypothesis: two.sided

testZeroInflation(mert.prop.simulation)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.98814, p-value = 1
## alternative hypothesis: two.sided

The Mertensia proportion data were not overdispersed or zero-inflated.

Next, we created the pollen limitation model for the Mertensia developed seeds. The model includes site and plot treatment as random effects like the models above. Now the proportion of developed seeds is the response variable.

mert$undev.seeds<-(mert$total.seeds-mert$dev.seed) #creating an undeveloped seed column for cbind

Mert.prop.treat.glmmtmb<-glmmTMB(cbind(dev.seed,undev.seeds) ~ treat + (1|site/plot.treat), family = binomial, data = mert) #model with pollen lim treatment and proportion of developed seeds using the glmmTMB package
summary(Mert.prop.treat.glmmtmb) #summary of model results

##  Family: binomial  ( logit )
## Formula:          
## cbind(dev.seed, undev.seeds) ~ treat + (1 | site/plot.treat)
## Data: mert
## 
##      AIC      BIC   logLik deviance df.resid 
##    770.7    780.8   -381.3    762.7       89 
## 
## Random effects:
## 
## Conditional model:
##  Groups          Name        Variance Std.Dev.
##  plot.treat:site (Intercept) 0.12831  0.3582  
##  site            (Intercept) 0.02192  0.1481  
## Number of obs: 93, groups:  plot.treat:site, 8; site, 4
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     0.62111    0.15662   3.966 7.32e-05 ***
## treatopen-hand -0.06824    0.07032  -0.970    0.332    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

We did not detect an effect of pollen supplementation on the proportions of developed seeds for Mertensia individuals.

We again checked to see if the data were overdispersed.

mert.prop.treat.simulation<-simulateResiduals(fittedModel = Mert.prop.treat.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
plot(mert.prop.treat.simulation) #plotting simulation

testDispersion(mert.prop.treat.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.0089, p-value = 0.952
## alternative hypothesis: two.sided

testZeroInflation(mert.prop.treat.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99128, p-value = 1
## alternative hypothesis: two.sided

Our pollen limitation model for Mertensia proportions was not overdispersed or zero-inflated.

Below are the plots for Mertensia proportion data.

mert.open<-mert.open[!is.na(mert.open$prop.dev.seeds),] #removing NA values from proportion column

ggplot(mert.open, aes(x=mert.open$relative.position, y=mert.open$prop.dev.seeds, group=plot.treat)) +
  labs(y="Proportion of Developed Seeds",x="Days Relative to Population Peak Bloom",title="Effect of Mertensia Blooming Time on Proportion of Developed Seeds") +
  geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
  scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
  geom_smooth()

ggplot(mert.open, aes(x=mert.open$deviation, y=mert.open$prop.dev.seeds, group= mert.open$early.late)) +
  labs(y="Proportion of Developed Seeds",x="Days Relative to Population Peak Bloom",title="Effect of Mertensia Deviation from Peak on Proportion of Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
  geom_smooth(method="glm",aes(color=early.late)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Before or After Peak") +
  scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))

Despite a marginally significant interaction effect, the effect of deviation on early vs. late individuals was not strong for Mertensia proportions of developed seeds.

ggplot(mert.open, aes(x=mert.open$number.conspecifics, y=mert.open$prop.dev.seeds, group= mert.open$plot.treat)) +
  labs(y="Proportion of Developed Seeds",x="Conspecific Density",title="Effect of Mertensia Conspecific Density on Proportion of Developed Seeds") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
  geom_smooth(method="glm",aes(color=plot.treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Plot Treatment") +
  scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))

Finally, we added pollinator activity.

ggplot(mert, aes(x=as.Date.factor(mert$midpoint), y=mert$prop.dev.seeds, group= mert$treat)) +
  labs(y="Proportion of Developed Seeds",x="Date", title="Pollinator Activity Period and Proportion of Developed Seeds in Mertensia") +
  theme_classic() +
  geom_point(size=3, alpha = 0.8, aes(color=treat)) +
  geom_smooth(aes(color=treat)) +
  scale_color_brewer(palette="Dark2") + 
  labs(color="Open or Hand-pollinated") +
  scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  annotate("rect",xmin=as.Date.character("2019-06-19"), xmax=as.Date.character("2019-07-03"),ymin=0,ymax=Inf,alpha=0.1) +
  scale_x_date(date_breaks = "3 days", date_labels = "%m/%d")

6.2 Delphinium

We had to find the total developed seeds and undeveloped seeds before creating the model. To do this, we combined the seed counts of the first flower in bloom with the rest of the flowers on the plant.

delph.open$dev.seed.total<-delph.open$dev.seed+delph.open$dev.seed1 #creating developed seeds column
delph.open$undev.seed.total<-delph.open$undev.seed+delph.open$undev.seed1 #creating undeveloped seeds column

#Delph.prop.glmmtmb<-glmmTMB(cbind(dev.seed.total,undev.seed.total) ~ deviation * early.late + plot.treat * number.conspecifics + rate + (1|site), family = binomial, data = delph.open) #modeling Delphinium proportion of developed seeds with glmmTMB package
#removed treat fixed effect (only open)
#summary(Delph.prop.glmmtmb)

Our model shows that relative position, plot treatment, and number of conspecifics had an effect on seed set. We see a marginally significant trend that late individuals have lower seed set than early individuals. The manipulated individuals had significantly higher seed set than the control individuals. Also, there was a significant interaction between number of conspecifics and plot treatment.

We again tested for overdispersion and zero inflation. We used the DHARMa package for the Delphinium proportion data.

#delph.prop.simulation<-simulateResiduals(fittedModel = Delph.prop.glmmtmb) #creating simulation for Delphinium model
#plot(delph.prop.simulation) #plotting simulation
#testDispersion(delph.prop.simulation) #checking for overdispersion
#testZeroInflation(delph.prop.simulation) #checking for zero inflation

Our Delphinium proportion data is overdispersed and zero-inflated. We adjusted our model to correct for these issues. To fix overdispersion, we are using a betabinomial family, and we added a zero-inflation formula to correct zero-inflation.

Delph.prop.glmmtmb<-glmmTMB(cbind(dev.seed.total,undev.seed.total) ~ deviation * early.late + plot.treat * number.conspecifics + rate + (1|site), ziformula=~1, family = betabinomial, data = delph.open) #modeling Delphinium proportion of developed seeds with betabinomial family and ziformula
summary(Delph.prop.glmmtmb) #summary of model output

##  Family: betabinomial  ( logit )
## Formula:          
## cbind(dev.seed.total, undev.seed.total) ~ deviation * early.late +  
##     plot.treat * number.conspecifics + rate + (1 | site)
## Zero inflation:                                           ~1
## Data: delph.open
## 
##      AIC      BIC   logLik deviance df.resid 
##    333.9    352.8   -156.0    311.9       30 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  site   (Intercept) 1.011e-09 3.179e-05
## Number of obs: 41, groups:  site, 4
## 
## Overdispersion parameter for betabinomial family (): 6.73 
## 
## Conditional model:
##                                       Estimate Std. Error z value Pr(>|z|)
## (Intercept)                         -5.550e-01  5.541e-01  -1.002    0.317
## deviation                            4.457e-02  8.726e-02   0.511    0.609
## early.latelate                       5.469e-01  7.269e-01   0.752    0.452
## plot.treatmanip                      5.551e-01  7.041e-01   0.788    0.431
## number.conspecifics                 -4.369e-05  8.929e-03  -0.005    0.996
## rate                                -1.014e-01  7.926e-02  -1.279    0.201
## deviation:early.latelate             8.399e-02  1.676e-01   0.501    0.616
## plot.treatmanip:number.conspecifics -2.342e-03  1.246e-02  -0.188    0.851
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -3.524      1.008  -3.494 0.000475 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

After correcting for overdispersion and zero-inflation, we did not detect an effect of relative blooming time, conspecific density, or plot treatment on the proportion of developed seeds for Delphinium individuals.

We again checked for overdispersion and zero-inflation in the new model.

delph.prop.simulation<-simulateResiduals(fittedModel = Delph.prop.glmmtmb) #creating simulation for Delphinium model
plot(delph.prop.simulation) #plotting simulation

testDispersion(delph.prop.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.1015, p-value = 0.376
## alternative hypothesis: two.sided

testZeroInflation(delph.prop.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.7123, p-value = 0.632
## alternative hypothesis: two.sided

The data are no longer overdispersed or zero-inflated.

We created the mixed effects model for the pollen supplemented Delphinium individuals. Like the Mertensia model above, site and plot treatment are random effects, pollen limitation treatment is the fixed effect, and proportion of developed seeds is the response variable.

delph$dev.seed.total<-delph$dev.seed+delph$dev.seed1 #creating developed seeds column
delph$undev.seed.total<-delph$undev.seed+delph$undev.seed1 #creating undeveloped seeds column

#Delph.prop.treat.glmmtmb<-glmmTMB(cbind(dev.seed.total,undev.seed.total) ~ treat + (1|site/plot.treat), family = binomial, data = delph) #modeling pollen lim and proportion of developed seeds with glmmTMB package
#summary(Delph.prop.treat.glmmtmb) #model summary

Delphinium proportions of developed seeds were significantly affected by pollen supplementation. Pollen supplementation significantly increased the proportion of developed seeds.

Below are our simulations for checking overdispersion and zero-inflation.

#delph.prop.treat.simulation<-simulateResiduals(fittedModel = Delph.prop.treat.glmmtmb) #creating simulation for Delphinium model
#plot(delph.prop.treat.simulation) #plotting simulation
#testDispersion(delph.prop.treat.simulation) #checking for overdispersion
#testZeroInflation(delph.prop.treat.simulation) #checking for zero inflation

The pollen limitation model for Delphinium proportion data is overdispersed and zero-inflated. As with the model above, we changed the family to beta binomial and included the zero-inflation formula.

Delph.prop.treat.glmmtmb<-glmmTMB(cbind(dev.seed.total,undev.seed.total) ~ treat + (1|site/plot.treat), ziformula=~1, family = betabinomial, data = delph) #modeling pollen lim and proportion of developed seeds with glmmTMB
summary(Delph.prop.treat.glmmtmb) #model summary

##  Family: betabinomial  ( logit )
## Formula:          
## cbind(dev.seed.total, undev.seed.total) ~ treat + (1 | site/plot.treat)
## Zero inflation:                                           ~1
## Data: delph
## 
##      AIC      BIC   logLik deviance df.resid 
##    654.0    668.6   -321.0    642.0       78 
## 
## Random effects:
## 
## Conditional model:
##  Groups          Name        Variance  Std.Dev. 
##  plot.treat:site (Intercept) 2.766e-08 0.0001663
##  site            (Intercept) 1.777e-01 0.4215895
## Number of obs: 84, groups:  plot.treat:site, 7; site, 4
## 
## Overdispersion parameter for betabinomial family (): 6.89 
## 
## Conditional model:
##                Estimate Std. Error z value Pr(>|z|)  
## (Intercept)      0.4443     0.2504   1.774    0.076 .
## treatopen-hand   0.1702     0.1764   0.965    0.335  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Zero-inflation model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -3.6522     0.7245  -5.041 4.63e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

After correcting the model, we did not detect an effect of pollen limitation on the proportion of developed seeds for Delphinium individuals. Again we checked for overdispersion and zero-inflation.

delph.prop.treat.simulation<-simulateResiduals(fittedModel = Delph.prop.treat.glmmtmb) #creating simulation for Delphinium model
plot(delph.prop.treat.simulation) #plotting simulation

testDispersion(delph.prop.treat.simulation) #checking for overdispersion

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted
##  vs. simulated
## 
## data:  simulationOutput
## ratioObsSim = 1.0861, p-value = 0.264
## alternative hypothesis: two.sided

testZeroInflation(delph.prop.treat.simulation) #checking for zero inflation

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.227, p-value = 0.848
## alternative hypothesis: two.sided